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Abstract 



We derive a formalism, the separation method, for the efficient and accurate cal- 
culation of two-body matrix elements for a Gaussian potential in the cylindrical 
harmonic-oscillator basis. This formalism is of critical importance for Hartree-Fock 
^ZJ • and Hartree-Fock-Bogoliubov calculations in deformed nuclei using realistic, finite- 

5 . range effective interactions between nucleons. The results given here are also relevant 

for microscopic many-body calculations in atomic and molecular physics, as the for- 
malism can be applied to other types of interactions beyond the Gaussian form. The 
" ^ . derivation is presented in great detail to emphasize the methodology, which relies on 

lO ! generating functions. The resulting analytical expressions for the Gaussian matrix 

^ ' elements are checked for speed and accuracy as a function of the number of oscillator 

shells and against direct numerical integration. 
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1 Introduction 



Gaussian interactions play an important role in the microscopic description of 
molecular and nuclear processes [T]. The Gaussian form represents a relatively 
simple tvi^o-body potential vi^ith a finite range, which is needed in many realistic 
descriptions of many-body systems. In nuclear physics for example, the Gogny 
interaction [2] 
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+^Wls ( Vi - V2) X 5 (fi - fa) ( V 1 - V2) ■ (ai + 0^2) 

+to (1 + xoP.) 5 (fi - rl) (^^^^^) + K^oui (1) 

where and P,- are spin- and isospin-exchange operators and p is the total 
nuclear density, gives the effective (in-medium) potential between nucleons. 
Two Gaussian terms appear explicitly with range parameters pi and 1x2- A 
spin-orbit term with strength Wls uses a Dirac-delta function, but extensions 
of the Gogny force have been proposed [3] that introduce a Gaussian form 
for this term. Finally the Coulomb interaction Vcoui ~ Vl^i "^2! between 
protons is clearly not of Gaussian form, but the mathematical framework 
presented in this paper can be applied equally well to a Coulomb potential. 

For the calculation of matrix elements in molecular, atomic, and nuclear 
physics, harmonic-oscillator functions provide a convenient and popular or- 
thogonal basis. The calculation of Gaussian matrix elements in a harmonic- 
oscillator basis, however, poses definite technical challenges in accuracy as 
well as execution time. In previous work [4], the separation method was in- 
troduced as a way of calculating the Gaussian matrix elements efficiently and 
accurately for systems with spherical symmetry. In the separation method, 
two-body matrix elements are expressed as a more manageable finite sum of 
products of one-body matrix elements. In this paper, we derive the separation 
method for a wider class of systems that exhibit axial symmetry. These results 
are crucial, for example, in microscopic calculations of nuclear fission using the 
Gogny force, where the nucleus elongates along a symmetry axis, until scission 
occurs. 

Fission calculations in particular bring to the fore many of the technical dif- 
ficulties involved in the computation of Gaussian matrix elements. On the 
other hand, microscopic calculations of fission using the interaction in Eq. ([T]) 
have had considerable success in recent years PfGfT] . and are therefore of great 
interest. In the microscopic description of fission, the matrix elements of the 
nucleon-nucleon interaction are typically used in a Hartree-Fock-Bogoliubov 
(HFB) procedure to construct a Slater-determinant wave function for the nu- 
cleus. Scission configurations are then found by driving the nucleus to such 
exotic shapes that the delicate balance between its surface tension and the 
Coulomb repulsion between the nascent fission fragments is broken. The proper 
identification of scission configurations and the calculation of their properties 
depend sensitively on accurate calculations of the matrix elements of the ef- 
fective interaction. Fission also implies the evolution of the nucleus through a 
variety of exotic shapes leading to scission. Therefore many sets of matrix el- 
ements need to be calculated, each set corresponding to a harmonic-oscillator 
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basis optimized for a particular nuclear shape, and each set requiring a large 
number of oscillator shells. The resulting large-scale computations can become 
very time-consuming and are prone to errors in accuracy. Thus microscopic 
fission calculations must rely on fast and accurate algorithms to evaluate the 
two-body matrix elements, such as the separation method. The separation 
method is especially well-suited to the HFB algorithm, because the coeffi- 
cients needed to calculate the two-body matrix elements derived in this paper 
can be calculated quickly once and for all, and stored with relatively little 
computer memory. 

The goal of this paper is to derive the separation-method formalism for Gaus- 
sian matrix elements in a cylindrical harmonic-oscillator basis, with particu- 
lar emphasis placed on the details of the derivation because of its relevance 
to other types of interactions, and other applications involving the harmonic- 
oscillator basis. In particular, we rely heavily on the power and versatility of 
generating functions to derive many of the present results. We also present 
the derivations in great detail because they are rather involved, and although 
the same results may be arrived at by alternate approaches, the formulas will 
tend to be much more cumbersome and less computationally efficient than 
the ones obtained by the generating-function methods outlined here. Because 
of the lengthy and detailed derivations involved, many of the intermediary 
results have been placed in the appendices. These intermediary results are 
important in their own right, as they provide useful properties of harmonic- 
oscillator functions in a cylindrical basis, and the mapping between cylindrical 
and Cartesian harmonic-oscillator bases. 

In section [2l the basic formalism for the calculation of both radial and ax- 
ial components of the Gaussian matrix elements by the separation method 
are derived. In section [3l the accuracy of the method is examined both rela- 
tive to direct numerical integration, and as a function of the number of shells 
in the oscillator basis. The execution times for the separation method are 
also compared to those of the numerical integration. The mapping between 
harmonic-oscillator function in polar and Cartesian coordinates, needed in 
the development of the separation-method formalism, is derived in appendix 
lAl In appendix iBl the Gaussian two-body potential, V {ri,r2), is written in 
separated form with respect to ri and r2- Formulas reducing the products of 
harmonic-oscillator functions are derived in appendix [Cl and provide a pow- 
erful tool in the evaluation of integrals involving those functions. In appendix 
[Dl the result quoted in [9] for the separation-method formalism in the case 
of large oscillator-shell numbers is derived in detail. Finally, in appendix [El 
we obtain a series expansion for the direct angular integral of the Gaussian 
potential, which we use in the numerical integration of the potential in section 

El 
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2 Theory 



2. 1 General formalism 



We wish to calculate matrix elements of the two-body potential function 



(2) 



in the cylindrical harmonic-oscillator basis. We will write the matrix elements 
as 

V,,u = {t3\V\kl) 

x1/ (fi, fa) {n; &±, h) $„(o_^(o^40 (^2; b±, h) (3) 

where we have introduced the stretched harmonic-oscillator basis functions in 
the cylindrical coordinates (p, zji^ 



^nr,A,n. {r; b±, h^) = $„,,,a (p, ip; b±) {z; 



nr,|A| 



(p; h 



(4) 



with the radial-component function 



defined in terms of associated Laguerre polynomials L\^j {rf) as a function of 
and with a normalization constant given by 



A/'„,,|A| = T- 



2nr\ 



.1/2 



K + |A|)! 



(6) 



The Cartesian, z-axis-component function in Eq. 



^ We will drop the qualifier "stretched" when referring to the deformed harmonic- 
oscillator function in subsequent discussion for the sake of brevity. 
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(7) 



is expressed in terms of Hermite polynomials Hn^ (^) with 



and normalization constant 



1 



The harmonic-oscillator functions defined in Eqs. ((41) and ([7]) satisfy the or- 
thonormalization conditions 



The parameters 6^ and hz appearing in the harmonic-oscillator function defi- 
nitions are usually treated as variational parameters in HFB calculations, and 
chosen to minimize the energy. 

The central idea in this paper is to express the two-body potential as a sum 
of products of one-body potential functions 



Then the two-body matrix elements can be written in terms of one-body 
matrix elements 



where we will show that this last sum is limited to a finite number of terms. 
It will be useful to separate the radial and Cartesian components in each 
one-body matrix element to write 







(8) 



5 



{i\fn,A,nJk) 



and, similarly, 



2lT 



pdp / c/(/? $ (,) (p, 6x) /„,,,A (p, 6j 

JO ' 



.1*', A(fe) 



(p, v'; &j 



X / rfz $ j>) (2; 6^) (z; b,) $ (2;; 6^) 
|/n,,A| A;) {i \fnjk) 



00 (•27r 

^ pdp dip ,^(,) (p, p; b±) ^nr,A (P, bj 



x<l> 



W.A(0 



71. 
00 



(p, v^; &j 



nr,A 



00 



so that we can write Eq. (I8|) as 



ijkl 



E(^l/n.,Al^) (j 

71r,A 









'^'n,,.,A 







(9) 



In the remainder of this section we calculate the explicit expressions needed 
to evaluate the matrix elements Vijki- 



2.2 Cartesian component 



(z) 

Here we derive an expression for the Cartesian component, V^^^.^ , in Eq. 
We will show that 



S = i/f^ E E ^,„«^),„a)/K,n.)(10) 



, (») (fc) 
mz = \n), —n). 



,2nz = \ni'' -niA,2 



where is defined by Eq. (IB.4I1 . the Z^j^^^ coefficients by Eq. (1C.6I1 . and the 
I{mz,nz) coefficients by Eq. ([TSll . 

We start by evaluating 
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dz $ 0) {z; h,) {z- h,) $ ro {z- b. 



Using Eqs. ( IB.lll which gives the explicit form of {z; b^) and Eq. (IC.ll l to 
reduce the product of harmonic-oscillator functions, 



^ r 



I \ ^ T-imz 



^2 <l>„^ (^;&^)$n, (2; b,) 



By orthogonality of the harmonic-oscillator functions this is simply 



fill 



where we must have 



ij) _ ^(0 < < n(^') + for the T^f,, coefficient 

^2 5^2 ^^^^^ 

to be non-zero. Next, we use the explicit form of fn^ {z; bz) from Eq. (|B.2l) to 
write 



[12) 



Two of the harmonic-oscillator functions can be replaced with a single one, 
thanks to Eq. flCH) . 



ybzVTT 



"z +"z 



^z )^z J— 00 



(i) (fc) 



The remaining integral, which we write in terms of the function 



(13) 



J —00 

where = Gy%zi can be calculated with the help of generating functions. 
Indeed, using Eq. (lA.ip to form the product of the harmonic-oscillator func- 
tions, we have for any ti and t2 
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m=on=o Vm!n! 

from which, multiplying by the Gaussian factors in the definition of / (m, n) 
and integrating both sides of the equation, 



*2 ^^g2il2/6,-22/b2+2i2z/B.~2VBz 
oo 

oo oo 2{™+")/2 



6,5,7r 5: 5: ^^t-t«J n) (14) 
m=on=o vm!n! 

The integral on the left-hand side can be evaluated by completing the square. 



^^^2tiz/h,-z^/hl+2t2z/B,-z^/Bl ^^t^/v 1'°° dze'^^'''^^^^ 



where we have defined 



_ 1 1 

"2 Z 

f— J- J z. 



Thus, the left-hand side of Eq. (fT4l l becomes 



-{b,t^-B,t2f/{yblBl) 



e 



which can be expanded as 



(-1)^+^ 

Comparing with the right-hand side of Eq. (fT4ll . we see that we must make 
the identifications m = 2p — q and n = g in order for the equation to hold for 
any ti and t2- Then, 
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LHS 



EE 

m=0 n=0 



2p 



\p+<i 



\q ) P 



m,2p—q 



and the comparison with the right-hand side of Eq. (fT4l) yields 



/ fm, n) 



^_^^(m+n)/2+n m + U 



Note that m + n must be even. We simplify this form further by noting that 



BIp = 1 + G, 
blu = l + Gz^ 



(15) 
(16) 
(17) 



where Gz is defined in Eq. ( IB. 41) . This leads us to write 



VI + V 2"*+" (^!22±ii)! (1 + G',)^"'+"^/^ 



m + n 



Some of the constants can be factored out by defining the coefficient 



/ fm, n) = 



/ (m, n) 



m\n\ 



, (m— n)/2 



(^)!(i + a 



.(m+n)/2 



m + n 



(18) 



Then, returning to Eq. ( fT3l ). we obtain after some simplification 



(19) 



Having derived the explicit forms in Eqs. (fTTl) and (fT9l) . we can express the 
Cartesian component in Eq. ([9]) as 
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I G,-1 



(i) , (fc) 



E 



,{i)_„(0 



where / {mz,nz) is given by Eq. (fTSll . and the T coefficients are given by Eq. 
( 1C.6I ). An alternate form of V^j^j^i was proposed by Egido et al. [9] which yields 
more accurate results for large oscillator shell numbers, and is derived as Eq. 
(tEH) in appendix [Dl 



2.3 Radial component 



A formula similar to Eq. ( fTOll can be derived for the radial component, V^J'jji, 
in Eq. We will show that 



r{r) _ ^± ^ rpn,- 

ij^l (7 , + 1 ^ ^ rS; 

1. „__n»i— n ^ 



G 



nn,-A(^)+A('=) 



nr=0 n=0 



' - A(') ;n''=' ,A(fe) nl-'' - AO) in''' ,A(0 



X/ [rir, -A(^') + A('); n, -A^ + A^^^) 



(20) 



where is defined by Eq. (IB. 101) . the T coefficients by Eq. (1C.9I) . and the 
/ coefficients by Eq. (l27l) . The indices nj; and n^ ^ are given by Eq. (IA.29I1 . 
where the bar indicates that — A*^-'^ and — A*^*) , respectively, should be used in 
that definition due to the complex conjugation in Eq. Q. 



Using Eqs. f lB.71) for the explicit form of $„,,,a (p, &±), Eq. (1C.7I) to reduce the 
product of harmonic-oscillator functions, and the orthogonality of harmonic- 
oscillator functions 



rir.A 



2tt 



pdp I dip $*o)^^(,) (p, b±) $„,,A (p, b±) 



n,-A(j)+^ 



1 "j.i 

2in,-A(j)+A(') 



X PC?P d(p $;_A(i)+A(0 (p, &±) $n,,A (p, &±) 

where the bar superscript in the nj i symbol serves as a reminder that we must 
use -A(J) in Eq. fTOOll . because of the complex conjugation. The condition 
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Snr<nj , comes about from the definition of the T coefficients in Eq. (1C.9I) . The 
other matrix element in the radial component of Eq. ([9]) is written explicitly 
using the explicit form for fnr,A {p, b±) in Eq. (IB. 81) as 



{i \fnr,A\ k) = pdp dip w_^(,) (p, ^; b±) /„,,a (p, ^; b±) 



and using Eq. (IC.7I) . the product of harmonic-oscillator functions can be re- 
duced 



I ■ \r I L\ _ -^-L-^2n^ + |A| „ _a{')+A('=) 

JO Jo 

x$n.,A (p, v^; -aw+ac^) (p, v^; 

where = G^^^fex, and the ^ in ^ is a reminder that we must use — A^*) in 
Eq. (IA.29I1 . The remaining integral to be calculated is 



I{n,M-n,M)^ r pdpT dipe-^'/i'^''-y^'/i''l) 
Jo Jo 

x'^ni.fci (p, B±) ^n2M (p, b±) (21) 

and can be evaluated using the generating function in Eq. (IA.3I) by writing, 
for arbitrary vectors ti and ^2, 



^-P^+2p-t,/B^-py{2Bl)^-tl+2p-t2/b^^py{2bl) 

/— oo oo 

fci=— oo ni=0 

/— oo oo 
Xn2,k2 (^2) '^'n2,fc2 (P^ b±) (22) 



note that, for clarity, we have explicitly written the parameter dependence for 
the normalization coefficients ^fnl,\kl\ {B±) and J^n2,\k2\ {b±) given by Eq. (I6|). 
Multiplying both sides of Eq. (l22l) by the Gaussian factor that appears in Eq. 
( I2T11 and integrating, we obtain on the left-hand side 
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oo i'2ir 



LHS = e-'^-'^ / pdp / d^e 
Jo Jo 

and on the right-hand side 



(23) 



oo oo oo oo 



RHS = -Blbl E E E E A/'n..|fc.|(5±)A/'n.,|fc.|(^J 

fci=— oo ni=0 A:2=— oo n2=0 

xx„i,fei (tl) Xn2,fe2 (^2) / (^1, h; n2, ki) 



(24) 



which contains the desired coefficients / (ni, fci; n2, ^2)- The integral in Eq. 
f l23ll can be evaluated by introducing 



V = 



B\ hi 
and completing the square, 



o , 9 9 /"CO /■27r 

LHS = e' '"-'^-'^ pdp rf(^exp 
z/ 

^^g(2fexBxri-r2-Bit?-biti)/(Bifei. 
V 



(^V^p- 



using Eq. (IA.27I1 with ti ti/ (fe^y^) and ^2 — ^ ^2/ iB±y/h'), this can be 
further expanded as 



n=0 fc=— 00 



Next, we use Eq. (|A.23P to eliminate the remaining exponential. 



LHS- 



2 1,2 00 00 00 00 



E E E E 



{n + mi)! (n + 1112)1 ,2 



=ofc=-oomi=om2=o mi\m2\ [n\] 



K\k\ {B±b, 



ti 



U 
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Using Eqs (IA.24I1 to eliminate the complex conjugation, and (1A.25I) to factor 
out the coefficients inside the x functions, this takes the form 

/ x-2(n+mi)-|fc| / N -2(n+m2)-|fc| 



XXn+mi-k (fl) Xn+m2,k (^2) 



Comparing this result for LHS with RHS in Eq. (l24ll for arbitrary vectors ti 
and ^2, we are led to conclude that 

/(ni,/ci;n2,fc2)=0 ifA:i + A;2^0 (25) 
We are also led to make the identifications 



n + nil = 
n + m2 = n2 
—k = ki 
k = k2 

which allow us to write 



^ 752 1,2 oo oo oo CXD rr, } 

LHS=^^Y. E E E ^ .KikiiB^b, 

2 „=ofc=-ooni=on2=o (^1 -^)' (^2 -n)! (n!) ' ^ 

-2ni-|fc| / \ -2n2-|fc| 



X 



(fe±v^) (5±v^) Xm~k (ii) Xn2,k (^2) (26) 



and therefore, assuming \ki\ = |/c2| = \k\ because of Eq. (l25l) . the comparison 
between LiJS* and RHS, in Eqs. fl26l) and f l24l) respectively, yields 



, 5,,+,„o(&xv^)-'"^-'"(i?xv^)"'"^-"^i!r.2! 
/ (ni, fci; n2, fc2) = 



■^m,\k\ iB±)J^n2,\k\ (bj 

00 I 

xE 



^0 (ni - n)! (^2 - n)\ (n!)^ 



xy^ni!(ni + |A;|)!n2!(n2+|/i;|)! 

00 1 

E = 

(ni - n)! (722 - n)!n! (n + \k\)\ 
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Using Eqs. (fT5ll - (|T7ll with G± instead of we can simplify the factor outside 
the summation 



B±b±iy 



-^2n2-\k\ 



(l + Gl' 



-ni-|fc|/2 



-n2-|fc|/2 



G 



gT + Gl'" 

(ni-n2)/2 



{gT + G 

and, for compactness of notation, we define 



1/2 l/2\"l+"2 + |A:|+l 



oo 2 

^ (m, ifci) = x: _ _ ^ 1^1^, 

which, after some simplification can be written as 



S(?ii,n2, \k\) 



^i!(^2 + |A;|)!^o 
1 

(ni +n2 + |A;|)! 



Til 

n 



( . 



n2 + \k\ 
\n J y n2 - n 

ni+n2 + \k\ \ I ni+n2 + \k\ 
rii \ n2 



where Eq. 0.156(1) in |8] was used to obtain the second line. Therefore, we 
finally have 



/ (ni, ki;n2, ^^2) = ^fci+fe.o 



G 



(n\-n2)/2 



^^1/2 _|_ ^-l/2^"l+"2 + |fc|+l 

x^ni\ {m + \k\)\n2\ {n2 + (m, ris, \k\) 



As in Eq. (|T8l) . it will be convenient to factor out some constant terms. There- 
fore we define 



I {ni,ki;n2,k2) = 



2ni + |fc| 



I (ni,ki;n2,k2 



^n,\in, + \k\)\n2\{n2 + \k\)\, 



and the radial component in Eq. Q becomes 



{n,,n2,\k\) (27) 



14 



nr. A 



OO OO "'i,k 



'-^-L nr=OA=-oon=0 

^nr,A'^'+A<') r r 
^ ^ ,A0) ;n« ,A(0 ^A,- AO) +A{0 

_ ± — ^ „„ _a{^)+A(*=) ,-aO)+a(o 

"(7, +1 2^ 2^ -AW;4'=',AW n^.-'' -AO);4'\a(0 

xl{nr, -A^^') + A^'); n, -A^*) + A^'^)) 

Thus , using Eqs. (fTOl) or f lD.2ll and (l20l) . the full matrix element V^jfci in Eq. 
([9]) can be calculated as an analytical expression. In the next section, we will 
examine the computational merits of these results. 



3 Discussion 



In this section, we will compare three different ways of evaluating the Cartesian 
i^ifki) radial {V^J'^i) components of the Gaussian matrix elements in Eq. 
([3]): 1) direct numerical integration of Eq. ((3]), 2) numerical evaluation of the 
separation-method equations (Eqs. ( fTOl l or ( ID. 21 ) for the Cartesian component, 
and Eq. (I20ll for the radial component) in double-precision mode, and 3) ex- 
act evaluation of the separation-method equations using the symbolic-algebra 
package Mathematica [10]. In principle, the first two methods-numerical eval- 
uation by either integration or the separation method-will give the values 
of V^i^j^i and V^J'j^i to within the limits of machine accuracy and roundoff er- 
rors, whereas the third-exact evaluation of the separation-method equations 
using Mathematica-will produce these matrix elements to any desired accu- 
racy (even beyond machine accuracy) and will serve as a reference check for 
numerical convergence of the integrals and roundoff errors. 

We begin by comparing the relative merits of the separation-method Eqs. 
( ITOll and (1D.2I) for the Cartesian component of the matrix element. The two 
equations are mathematically equivalent, but Eq. f lD.2l) was obtained from Eq. 
(fTOl) specifically to provide greater accuracy in numerical calculations. For all 
quantitative applications in this work, we have used 

/i = 1.2 fm 
bz = 3.3 fm 
6± = 2 fm 

These values of n, bz, b± are typical in HFB calculations using the Gogny 
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interaction for ^"^^Pu along the most likely path to scission [TT] . 



In practice, both Eqs. (fTOll and (ID.2I1 can be evaluated efficiently because the 
T^^ „2 ^iid I (m, n) or F^_^ coefficients can easily be calculated once and for 
all and stored with relatively little memory, to be used in reconstructing the 
matrix elements V^^^^ whenever they are needed. However, for large values 
of the quantum numbers n^, n^, n/c, and n; the sums in Eq. (fTOl) rapidly 
lead to sizable numerical inaccuracies. These inaccuracies arise because the T 
coefficients grow progressively larger with increasing values of the arguments, 
whereas the / coefficients decrease. The resulting sum of products of small 
and large numbers in Eq. (fTOl) becomes numerically unstable. The formula 
obtained by Egido et al. in [^, and derived as Eq. (ID.2P in the present work, 
avoids this problem. 

Fig. [H gives the maximum deviation between matrix elements calculated using 
numerical evaluations of Eqs. (fTOl) and (ID. 21) . To generate the plot, the equa- 
tions were compared for calculations of V^^^; as a function of the maximum 
harmonic-oscillator shell number A^o? i-C- for all possible quantum numbers 
such that < ni,nj,nk,ni < Nq, and the largest deviation was recorded for 
each point on the plot. We will refer to Nq as the size of the basis in the dis- 
cussion below. The deviations plotted in Fig. [Dare based on the dimensionless 
Gaussian function in Eq. ([2]), but with realistic interaction strengths for the 
Gogny force [l2], a deviation as small as 10~^ on the plot, can correspond 
to a discrepancy of the order of an MeV. Thus, for Nq greater than about 
16, Eq. (ID. 21) should certainly always be used instead of Eq. ( ITOl l. and in the 
remainder of this paper we will use it consistently for all A^o instead of Eq. 

Next, we compare an exact evaluation of Eq. (|D.2l) to the numerical inte- 
gration of the Cartesian component in Eq. ([3]). We choose to compare the 
separation method to a numerical integral of the potential because the latter 
is easily implemented, requires very little computer memory, and can be made 
arbitrarily accurate. The exact evaluation of Eq. ( ID. 211 was obtained using the 
symbolic-algebra package Mathematica. Within Mathematica, the expression 
in Eq. (ID. 21) was first reduced by symbolic manipulation to the exact algebraic 
form aVb/c, where a, b, and c are integers, for each choice of the quantum 
numbers rij, rij, Uk, and ni. That algebraic number could then be evaluated 
numerically to any desired accuracy. The numerical integration, on the other 
hand, was performed by Gauss-Hermite quadrature in double-precision mode 
(i.e., with 16 significant figures). The purpose of the comparison between the 
exact evaluation of Eq. ( ID. 21 ) and the numerical integration is to show that 
the numerical integration can be made arbitrarily close (up to the limits of 
machine accuracy) to the exact result, thereby validating Eq. ( ID.2p . In Fig. 
[2l the maximum deviation between the exact calculation and numerical in- 

(z) 

tegration of the V^jj^i values is plotted as a function of the number Nq^^^ of 
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iz) 

Figure 1. Maximum deviation between calculations of the matrix elements V^^^^ using 
the separation method in Eq. (fTOl) on one hand, and Eq. (ID.2|) on the other, plotted 
as a function of basis size A'o . 

quadrature points for a basis size A^o = 12. For N^faaA > 208, the limits of 
machine accuracy are reached in the numerical integration, and the maximum 
deviation betv^^een the tv^^o methods of calculating vj^^^i matrix elements levels 
out slightly above 4.3 x 10"^^. 

In Fig. m we compare the exact evaluation of Eq. (ID. 21) using Mathematica to 
its numerical evaluation in double-precision mode, as a function of basis size 
Nq. The trend in Fig. [3] shows the effect of roundoff error in the numerical 
evaluation of Eq. (ID. 21) . However, despite a clear decrease in accuracy with 
increasing basis size. Fig. [3] shows that a double-precision numerical evaluation 
of Eq. (ID. 21) still gives the value of the matrix elements vj^^^i to a very high level 
of accuracy. Even for a basis size as large as Nq = 24,, the largest deviation from 
the exact values is still only 1.5 x 10~^. For the remainder of this discussion, we 
will use the numerical evaluation of Eq. (ID.2I1 in double-precision mode rather 



than the exact Mathematica result, because the Mathematica calculations are 
prohibitively time-consuming, and the accuracy of the numerical evaluation of 
the separation-method formulas is more than sufficient for most applications. 

In Fig. m we extract the number of Gauss-Hermite quadrature points required 
by the numerical integration to obtain values that are satisfactorily close (say, 
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Figure 2. Maximum deviation between the numerical integration of the matrix el- 
ements V^^^^i and their exact evaluation using the separation method in Eq. (|D.2p 
with Mathematica for basis size A'^o = 12, plotted as a function of the number of 
Gauss-Hermite quadrature points in the integral. 

within a 10~^ discrepancy at most) to the values given by a numerical eval- 
uation of Eq. f lD.2[ ). The number of quadrature points plotted as a function 
of basis size A^o is moderately large, and increases steadily vi^ith Nq. Further 
below we will gauge the cost in computational time incurred by the numerical 
integration with these relatively large numbers of quadrature points. 

(r) 

We carry out a similar analysis for the radial component, 1^^^^;.;, of the matrix 
elements. In this case, for a given basis size Nq, the quantum numbers for 
the radial matrix element v}j^i in Eq. ([3]) take on all values such that < 
2nr + |A| < A^o with > 0. As we did in Fig. [2]for the Cartesian component, 
we compare in Fig. [5] an exact (Mathematica) calculation of Eq. (l20l) to a 
numerical integration of the radial component in Eq. ([3]) using double-precision 
Gauss-Laguerre quadrature, for a basis size A^o = 8. In Fig. El the maximum 
deviation between exact evaluation and numerical integration, plotted as a 
function of the number A'^quad of quadrature points, is made arbitrarily small 
with increasing N^y^^^ values until the limits of machine accuracy and roundoff 
error are reached for A^quad ^ 48, where the maximum discrepancy settles 
above 1.3 x 10"^^ 



18 




Figure 3. Maximum deviation between the numerical calculation and exact Mathe- 

(z) . 

matica evaluation of the matrix elements V-jj^^ using the separation method in Eq. 
(lD?2ll . plotted as a function of basis size A'^o- 

A comparison between exact (Mathematica) and double-precision numerical 
evaluations of the separation-method result in Eq. (i20l ) is plotted in Fig. [6] as 
a function of basis size Nq. The accuracy of the numerical evaluation clearly 
deteriorates with increasing basis size, but remains quite good nevertheless, 
reaching only a 1.2 x 10~^ maximum deviation for Nq = 12. For practical rea- 
sons, we will use the numerical evaluation of Eq. (l20l) in the remainder of this 
discussion, rather than the exact-but much slower-Mathematica calculation. 

The number of Gauss-Laguerre quadrature points needed to obtain a discrep- 
ancy of 10~^ or less between the numerical integration and numerical separa- 
tion method for v}j',ji matrix elements is plotted in Fig. [7] as a function of basis 
size. As in Fig. [4] for the Cartesian matrix elements, the required number of 
quadrature points is moderate and increases with basis size. The impact of 
these numbers of quadrature points on execution time will be investigated 
next. 

We now compare execution times for the numerical integration and numerical 
separation methods. The numerical integrations for the Cartesian and radial 
components are performed with the number of quadrature points given in 
Figs. [4] and Ui respectively, to ensure agreement to 10~^ or better with the 
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Figure 4. Minimum number of Gauss-Hermite quadrature points needed to achieve 
10~^ or better agreement between the numerical integration of the matrix elements 
V^^-^^; and their evaluation using the separation method in Eq. (ID.2p . plotted as a 
function of basis size A'^o . 

separation-method results. In order to speed up the numerical integrations, 
the harmonic-oscillator functions are calculated at the appropriate quadrature 
points and stored once and for all. A set of nested loops then evaluate the 
multidimensional integrals by recalling the stored values of the functions as 
the terms in the quadrature are summed. Likevi^ise, for the calculations by the 
separation method, the T, /, and F coefficients are calculated ahead of time 
and recalled as needed in the evaluation of the matrix elements using Eqs. 
dnj) and ((201). 

The calculations have been performed on a 2.13-GHz Pentium M proces- 
sor in double-precision mode. The execution times are plotted in Fig. [8] for 
the z component of the matrix element, and in Fig. [9] for the radial compo- 
nent. The times plotted include the setup time needed to pre-calculate the 
harmonic-oscillator function values and separation coefficients appropriate to 
each method. The difference in execution times between the numerical and sep- 
aration methods become staggering with increasing basis size. For large-scale 
computations requiring matrix-element calculations over a range of values of 
the harmonic-oscillator parameters and 62, such as maps of fission shapes 
for a single nucleus or maps of nuclear properties for large sets of nuclei, direct 
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Figure 5. Maximum deviation between the numerical integration of the matrix el- 
ements V^^^^ and their exact evaluation using the separation method in Eq. (I20|l 
with Mathematica for basis size A'^o = 8, plotted as a function of the number of 
Gauss-Laguerre quadrature points in the integral. 

numerical integrations rapidly become unfeasible without parallel machines. 
Even with parallel processing, modern nuclear-physics problems (e.g., the 
microscopic treatment of fission in a multidimensional collective-coordinate 
space) will eventually overwhelm any given computational resource, and in 
order to match the accuracy of the separation method, numerical integrals 
will generally require an inordinate number of quadrature points. 



4 Conclusion 

We have derived explicit expressions for Gaussian matrix elements in a cylin- 
drical harmonic-oscillator basis, using the separation method. These expres- 
sions have been tested against direct numerical integration and found to be 
highly accurate and computationally efficient. These characteristics make the 
separation method an invaluable tool for computationally-intensive applica- 
tions, such as the microscopic description of fission. The work presented here 
has wider relevance than to the Gaussian form, or to nuclear-physics problems 
alone. In particular, the methodology used in the present derivations, which 
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Figure 6. Maximum deviation between the numerical calculation and exact Mathe- 

(r) . 

matica evaluation of the matrix elements V-jj^^ using the separation method in Eq. 
(I2?l . plotted function of basis size A'^o . 



relies heavily on generating functions, can be applied to other types of inter- 
actions and a wider class of basis states to derive analytical, computationally- 
efficient expressions for matrix elements. For example, in future publications, 
we will apply the separation method to the Coulomb and Yukawa interactions, 
and extend the formalism to bases of displaced and two-center deformed har- 
monic oscillators. These planned extensions to the separation formalism en- 
large the range of applications of the method to many problems of central 
importance in nuclear, atomic, and molecular systems. 



We wish to thank D. Gogny for invaluable guidance in the development of 
the formalism and preparation of this manuscript. This work was performed 
under the auspices of the U.S. Department of Energy by Lawrence Livermore 
National Laboratory under Contract DE-AC52-07NA27344. 
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Figure 7. Minimum number of Gauss-Laguerre quadrature points needed to achieve 
10~^ or better agreement between the numerical integration of the matrix elements 
V-jf^l and their evaluation using the separation method in Eq. (1201) . plotted as a 
function of basis size A'^o . 



A Mapping between Cartesian and polar coordinates for harmonic- 
oscillator functions 



In this section, we derive an identity relating the harmonic-oscillator functions 
expressed in two-dimensional Cartesian coordinates (x, y) to those in polar 
coordinates (p, ip) where 



2 2,2 

p =x + y 
y 

tan ip = — 

X 

To this end, we will first need to derive generating functions for the harmonic- 
oscillator functions in the two coordinate systems. 
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Figure 8. Comparison of total execution times for the evaluation of F^^^ by numerical 
integration and by the separation method in Eq. (ID. 21) . as a function of basis size 
iVo. 

A.l Generating function in Cartesian coordinates 



In this appendix, we derive the generating function 




(A.l) 



for the Cartesian harmonic-oscillator functions in Eq. (JT]). 



We begin with the generating function for Hermite polynomials (Eq. 8.957(1) 
p. 1034 in [8]), for arbitrary variables x and t, 



oo j.k 

making the substitution x ^ x/hm order to introduce the harmonic-oscillator 
parameter 6, 
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Figure 9. Comparison of total execution times for the evaluation of F^^^ by numerical 
integration and by the separation method in Eq. (f20]l . as a function of basis size A^o- 



-t^+2tx/b . 



Next, we introduce the Gaussian and normalization factors appearing in the 
definition of the harmonic oscillator function in Eq. ([7]) 



oo 



or, in terms of the harmonic-oscillator functions, 
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A. 2 Generating function in polar coordinates 



Here, we derive a generating function for the polar harmonic-oscillator func- 
tions defined in Eq. dH), 



, , oo oo { _T +2n+\k\ 

^ o.Jn\ {n+\k\)\ 



(p, ip; h) 



k=—oo n- 



(A.2) 



which we also cast in the form 



'P+2p-t/b~p'^/(2b^ 



/— oo oo 
^ I , n 



fc = — oo n=0 



(A.3) 



where the functions Xn,k (tj are defined by Eq. (IA.9I) . 

To derive a generating function for harmonic-oscillator functions in polar coor- 
dinates, we begin with the generating function for Laguerre polynomials (Eq. 
8.975(3), p. 1038 in [Hj), for arbitrary variables x and and a > — 1 



oo 



„t^or(n + a + l)^"^'^^ 



(A.4) 



In order to match the definition of the harmonic-oscillator function in Eq. (I5|), 
we substitute \/x = p/b, ^Jz = —it, and a = \k\ where k is an integer. Then, 
isolating the Bessel function on the left-hand side, Eq. (1A.4I ) takes the form 



.,K-2>M/.)^e.^-0l''(f)';||^.-(g) (A.5, 

On the other hand, the generating function for a Bessel function of the first 
kind for arbitrary z and (p is (Eq. 8.511(4), p. 973 in [8]) 



^iz cos j2 i\jj^{^z)e'''^ 

fc=— oo 
oo 

= ^'"^|fc|(^)e^'^ (A.6) 

fc=— oo 

where the second line follows from Eq. 8.404(2) in [8]. Substituting z = 
-2ipt/b'^ into Eq. (IATbI) . 
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g2ptcos^/6_ j2 (_2i^) e'^^ (A.7) 

fc=— oo 



b 

Finally, plugging Eq. ( 1A.5I1 into Eq. ( 1A.7I) yields 



fc=— oo \ b J 



where the right-hand side can be made to look more like the harmonic- 
oscillator function definition in Eq. (I5|), 



oo oo 



e 



t +2pt cosip/b . 



E E(-i) 



^2n+|fc| y2^e^'/(''') 



fc— ^— {n+\k\)\ Afn,lk\ 

X 



,6/ " VV VSF 

or, after straightforward simplifications. 



^-*^+2,tcos,/.-pV(2.^) = 6v/^ E E $n.. (P, b) 

fc=-oon=o yn! (n + \k\)\ 

Note that there is a potential ambiguity in the meaning of the angle (p in 
Eq. ( IA.2II . In fact, Eq. (IA.2I 1 was derived for any arbitrary value of ip but on 
left-hand side, the term pt cos p in the exponent suggests a dot product p ■ t 
with (y9 the angle between the vectors, while on the right-hand side, writing the 
harmonic-oscillator function ^ (p, V^; b) suggests that p is the polar angle of 
the vector p. To lift this apparent ambiguity, we introduce the polar angle (ft 
of vector t explicitly by noting that if 9 is the angle between vectors p and t 
with 6 = if — ift, then according to Eq. (I4|) 



$n,fc (p,^;6) = <l>„,|fc| (p; &)^= 
and therefore 

(p, 6) = <l>„,fc (p, 6) e-^'^'^* (A.8) 
Writing the left-hand side of Eq. ( 1A.2I1 in vector form, we now have 
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, . oo oo "l" +2n+|fc| 

fc=-oon=o yn! (n + |A;|)! 



For convenience, we introduce the function 



X„,,(t) = ^^t2"+l'=le-^'=^' (A.9) 

which allows us to write the generating function for polar harmonic-oscillator 
functions as 



' ^ fc=-oon=0 

This form will be convenient for some derivations, and we will obtain useful 
properties of the function Xn,k in section lA.Si 



A. 3 Polar-to- Cartesian mapping 



Having derived generating functions for the harmonic-oscillator functions in 
both polar and Cartesian coordinates, we can now obtain a relation between 
the two, 



rix+Uy 

<|.„^ (x; b) {y- b)= E C:T'^ ^.^r.,-m . (P, y^; b) 



(A.IO) 



where the coefficients C"^^'"^^ are given by Eq. (lA.lTp . 

In order to relate the polar and Cartesian harmonic-oscillator functions we will 
use Eqs. (lA.ip and (lA.2p . We will assume axial symmetry and use the same 
parameter h for all the coordinates involved. Consider the arbitrary vectors 
p = XX + yy and t = t^x + tyj) in the two-dimensional Cartesian coordinate 
system, with p-t = pt cos 9. Note that we are using the symbol 9 for the angle 
between vectors p and t. We can write 



^-tl+2t^x/b-x^/(2b'^)^-tl+2tyy-y'^/(2b'^) ^ ^-t^ +2pt cos e/b- p'^ / (2b'^) 

Using Eqs. flA.ll) and (IA.2I ). this can also be written as 
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oo oo ryinx+ny)!"! 



oo 



°o ( — ]]^f'^n+\k\ 



k=—oo n- 



1o^n\ {n+\k\)\ 



(A.ll) 



We must now equate the terms on the left-hand side to those on the right- 
hand side. We would like to introduce the polar angle ip of the vector p instead 
of the angle 9 between p and t in these expressions, because the final result 
should be completely independent of the choice of vector t. Using Eq. ( 1A.8I) . 
Eq. (lA.llI) becomes 



oo oo n{nx+ny)/2 /j- \nx / + 



(x; b) {y- b) 



oo oo 

^ E E 

fc=— oo n 



2n+|fc| 

^n! {n+\k\y. \bJ 



(A.12) 



All we have to do now is identify terms on the left- and right-hand sides. We 
can establish this correspondence by expressing t and ipt in terms of and ty. 
To this end, we write 



t 



(te" 



\k\ 



where we have introduced the sign quantity 



Sk = 



-1 



A; > 
A; < 



(A.13) 



Note that we can write 

^^-iskft _ ^ _ it gin (sk^pt) 

= t cos (pt — iskt sin ipt 

where the second line follows because Sfc = ±1. Thus we have 



^2n+\k\^-ikipt 



2 2\ / 



^EE 

p=0 13=0 



n 



I 



k\ 



V 1 
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We substitute this result into the right-hand side of Eq. (IA.12P to get 



oo oo 



RHs= E E 



fc=-oon=o ^Jn\ {n + 




-ISk 



,\k\-q 



(A.14) 



Comparing with the left-hand side of Eq. ( 1A.12I ). we see that we will need to 
make the identifications 



2p + q = n^ 
2n + \k\ — 2p — q = ny 

which also implies the important relation 



+ n„ = 2n + I A; I 



(A.15) 



We wish to replace the sums in Eq. (IA.14P over n and p with sums over 
and Uy. Since = ^p + q, it is clear that will span the full range of integers 



starting with 0. Similarly, Eq. (IA.15I) implies that 



2n + \k\ — and for 



any n^j there will always be a set of n and k values such that Uy spans the 
full range of integers from 0, independently of the value of index rix. Thus we 
can make the substitution 



oo oo 



EE-E E 

n=Op=0 nx=Ony=0 



Next, we note that Eq. ( 1A.15I 1 can also be written as 2n = + Uy — \k\, and 
since n > 0, we must therefore have \k\ < + Uy. Finally, 2p = — q, and 
since p > 0, we conclude that q < n^. Thus we can also make the substitution 



oo \k\ rix+Uy min(n^,|A;|) 

E E- E E 

fe=— oo q=0 k=—nx~ny q=0 



and Eq. (1A.14I1 becomes 



rix+Uy rmn{nx,\k\) 

E E 



Jnx+ny-\k\)/2 




nx+ny-\k\ I nx+ny + \k\ \ 
2 ■ 2 ■ 



2 



(A.16) 
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Note that in the sum over k, the index can be stepped by 2 units at a time, 
because of the restrictions imposed by the factorials. Comparing the left-hand 
side of Eq. ( 1A.12I ). and its right-hand side given by Eq. f lA.161) . we deduce 



Tlx • • 



(x; b) <i>n, [y] h) 



-isk) 



rix+riy min(na;,|fc|) / -, \(nx+ny-\k\)/2 

y y L^i 

/ nx+ny-\k\ 



\k\-q 



\ 



Hx-q 

2 



\k\ 



*^*nx+ny — |fc| ^ 



or, in more compact notation. 



(x; b) <!>ny [y] h) : 



rix+Uy 

^ _ E 



where 



'^n,k 



njnyl (_i) 



{na:+ny- |fc|)/2 min(rai,,|fc|) 



2{na:+ny)/2 l nx+ny-\k\ ^ na:+ny+\k] ] ^ 
V 2 ■ 2 ■ 



V 



2 

rix — q 



\k\ 



(A.17) 

The appearance of the index n in the symbol C^^fc"^, even though it is not 
explicitly used, serves as a reminder of the implicit relation between the indices 
given by Eq. (1A.15I1 . 



A. 4 Cartesian-to-polar mapping 



In this section, we derive the inverse transformation corresponding to Eq. 

(lAloD, 



2n+|fc| 

'^n,k (P, V;b)= E Clll^,ny'^2n+\k\-ny {x; h) <l>„^ (?/; h) 



(A.18) 



which expresses the polar harmonic-oscillator functions in terms of the Carte- 
sian functions. The coefficients C^'^^^ are given by Eq. (IA.2ip . 



We start again from Eq. (IA.12p . but this time, we express t^ and ty on the 
left-hand side in terms of t and ipt. Consider then 
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+ e 



-upt 



2i 



Expanding the powers and grouping terms yields 



-j-Ux+riy Tlx 

tyty = _ — _ 

p=0 q=0 



-x y 



Substituting into the left-hand side of Eq. (IA.12P produces 



oo oo 



LHS= E E 



Tlx 

EE 



^^-iinx+ny-2p-2q)^,^^^ (x; h) {y, h) 



\ p 



■^riy-q 



Comparing with the right-hand side of Eq. (IA.12I1 we see that we need to make 
the identifications 



n^ + ny = 2n+\k\ (A. 19) 

nx + ny-2p-2q = k (A. 20) 

we therefore introduce a summation over n and k with the help of Kronecker- 
delta functions, 



oo oo 



LHS=Y E t2"+|fcle-^'='/'t2- 



{nx+ny)/2 



E E^^ 



nx+ny,2n+\k\ 



n=0 fc=— oo 



X . 2^ 2^ 02p+2q,nx+ny~k 



n^\ny\i'-y p=oq=o 



(-1)"^ 



Uy-q 



where the Kronecker-delta functions collect those terms in the remaining sum- 
mations needed to satisfy Eqs. (lA.lOp and (lA.20p . The restrictions imposed by 
the Kronecker-delta functions can be used to eliminate the summations over 
Ux and p 
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oo oo 



l^Ug — ^ ^ ^2n+\k\^-ikipt2-n-\k\/2 ^ ' ^2n+\k\-ny {^'i b) ^riy {U'l b) 



n=Ok= 



min{ny,n—q+(\k\ — k)/2) 

E 

9=0 



2n + \k\ — n 



\n-q+^ J yq 



^0 {2n + \k\ - ny)\nyW 



Comparing with the right-hand side of Eq. ( 1A.12I1 we deduce the relation 



2-n-|fe|/2 ^'y^^ *2„+|fc|-„, (x; b) {v, b) 



ny=o \J{2n + \k\ - nyy.Uyli'^y 



Qmax 

q=0 



2n + \k\ — n 



1 


( ny\ 




\ 9 / 



-1 



^/ri! (n+ 



where 



gmax = niin (rij^, n + {\k\ - k) /2) 
which we write as 



^n,k (p, v^; &) 



2n+|A:| 



(p, 6) = ^ (:7;^;^„^$2„+|fc|_„^ (x; b) <l>„^ (y; 6) 



r!,a=0 



with 



/^n,k 
nx,ny 



1 









\ny-q 



(A.21) 

The appearance of the index Ux in the symbol C^^^„^, even though it is not 
explicitly used, serves as a reminder of the implicit relation between the indices 
given by Eq. (IA.15I) . 



A. 5 Properties of the function Xn,k (tj 



In section IA.2I we introduced the function Xn,k (tj which was used to obtain 
a generating function for harmonic-oscillator functions in polar coordinates. 
This function has many useful properties which we will exploit in further 
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derivations. In this section, we obtain some important properties of Xn,k [t 
From the definition of the Xn,k (t) function in Eq. (KM), 



^ ^ nl 
we can easily show that 




As a corollary, we can use Eq. ( 1X221 ) to show 



m=0 



m\n\ 



The complex conjugate of Xn,k is also readily expressed as 



(A.22) 



(A.23) 



Xl,k {t) = Xn,-k (t) 



(A.24) 



and a scale factor can be factored out. 



n,k {o,tj 



(A.25) 



Next, We will use the function Xn,fc, to expand the expression exp {2ti ■ t2^. 
Starting with the generating function for Bessel functions of the first kind, Eq. 
( 1A.6I1 with z = —2itit2 and (p = ipi ~ ip2, 



(A.26) 



k=—oo 



Next, we use the series expansion for Bessel functions (Eq. 8.440 in [8]), 



^.W-^J ^^k\{y^k)\ V2 
to write Eq. ( 1A.26[ 1 as 
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k=—oo 



n=0 



nl fc + n ! 



oo oo 



n?o.Ho^!(|A:|+n)! 



(tits 



\2ri+|A;| ik{ifii-ip2) 



or, 



1.2 oo oo 
n=0 k=—oo 



^l\k\ iP) X*n,k ifl) Xn,k {h 



(A.27) 



where ^fnr,\A\ (b) is given by Eq. ([6l), and the oscillator parameter b cancels out 
in the right-hand side. Next, we derive an expression for the product of two 
Xn,k functions, using the definition in Eq. f lA.911 



Xni,ki 



-1) 



ni+n2 



_^2ni +2712+1 fci| + |fc2 1 g-«(fcl+fc2)¥'t 



ni\n2\ 

at this point, it is convenient to define the quantities 



(A.28) 



nio = ni +n2 + 



\ki\ + l/cal - \ki + k2\ 



\ki\ + \k2\ - \ki + k2\ 



(A.29) 
(A.30) 



which recur throughout the paper. Then Eq. (1A.28I) becomes 



Xni,fci (t) Xn2,fc2 (^) = (-1) 



-fci,2 ^1,2! (-1) ''" ^2ni.2 + |fci 

ni\n2\ nio\ 



^zni,2-|-|fi;i+/c2|g-j(fci+fc2)</3t 



or, 



(A.31) 



Next, we obtain an expression for the function Xn,k {ti + ^2) of a sum of vec- 
tors. We write for an arbitrary vector t 



g2(ti+t2j-t ^g2ti-tg2t2-t 



(A.32) 



Using Eq. ( 1A.27I1 . the left-hand side is 



r 2 CO 00 

LHS=-Y. E Km (b) Xl,k {ti + ^^2) Xn,k (t) 



(A.33) 



n=0 fc=— 00 
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while the right-hand side of Eq. (1A.32I) is 



l4 oo oo oo oo 



ni=0 fci=— oo n2=0 fc2=— oo 



Using Eq. (lA.Sip . this reduces to 



rA oo oo oo oo 

RHs=-Y. E E E ACVi(^X,i^.i(^)(-i) 



fcl2 ^1,2! 



ni=0 fci=— 00 n2=0 fc2=— 00 



?T,i!n2! 



In order to compare with Eq. f lA.331 1. we introduce summations over the indices 
n and k with the help of Kronecker-delta functions, 



t4 00 00 00 00 00 00 

RHs=-Y. K,\^.\iP)ui,,^,{h){-r 

n=0 fc=— 00 ni=0 fci=— 00 n2=0 /c2=— 00 
^Xni.fci (^1) Xn2,fc2 (^2) ^n,ni,2^kM+k2Xn,k {t) 



ki,2 '^1,2! 



ni\n2\ 
(A.34) 



Comparing Eqs. (IA.33P and (IA.34I1 for an arbitrary vector t, and taking the 
complex conjugate, we are lead to write 



n,k [tl + ^2) — E 



00 00 00 00 



E 

ni=0 fci=— 00 n2=0 ^2=— 00 



E E E 



DZtkv,n2M^'^^M i^l) Xn2,k2 (^2) 



where 



D 



n,k 

ni,fci;7i2,fc2 



-, \ni+n2"n ip ~^ 1^1)- r r 

{ni + \ki\)\ {n2 + \k2\)\ 



(A.35) 



(A.36) 



Note that we have used the condition imposed by the Kronecker-delta function 
5n,ni,2 ^iid the definition of ni 2 in Eq. (|A.29l) to write 



:-i)'^"=(-i) 



Tn+n2—n 



Finally, we derive an expansion for the product exp {2ti ■ ^2) Xn,k {ti + ^2) • 
Though it is tempting to use Eq. (IA.27I1 for this, we will adopt a different 
approach which will yield a simpler expression in the end. We write 
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We treat the first exponential on the right-hand side using Eq. ( 1X231 ). so that 



Next, we use Eq. (IA.35P to expand the Xn+m,k \ ti + ^2) function 



[-ir{n + m)\ 



^ ' m!n! 



m=0 

00 00 00 00 

X V V V V n'^+'"-'= 

Z.^ Z.^ Z.^ Z.^ -^ni,fci;n2,A:2 
ni=0 fci=— 00 »i2=0 A;2=— 00 

XXni,fci (^l) Xn2,fc2 (^2) 



and use Eq. (IA.23I1 again to eliminate the remaining exponential on the right- 
hand side 



^ ^ 00 00 00 00 00 00 

e^*''*^Xn,fc (?L + ^2) = XI XI XI XI XI XI -^nl^fci,mi;n2,fe2,m2 
mi=0 m2=0 ni=0 fci=— oo n2=0 fe2=— oo 

^Xni+mi,fei (''l) Xn2+m2,k2 {h) (A.37) 

where we have defined 

^ >^ (-1)" {n + m)\ {ni + mi)! (ng + ma)! n+m,fc 

ni,fci,mi;n2,fc2,m2- ni,fci;n2,fc2 

m!n!mi!ni!m2!n2! 

{n + m+\k\)\ ^ 
(r^i + |fci|)!(n2 + |fc2|)! '''^^'^ 

which simplifies to 

k _ ^^1,2! (^1 + "^1)! (^2 + ^2)! (ni,2 + + A;2|)! 



D 

ni,ki,mi;n2,k2,m2 



(ni,2 - n)!n!mi!ni!m2!n2! {rii + |A;i|)! (n2 + |A;2|)! 

X ( — l)"i+"2 " 5„<„j 2^fc,A:i+fc2 (A. 38) 



Note the disappearance of the infinite sum over m in favor of the Kronecker- 
delta function 6n<ni 2 ■ 
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B Decomposition of two-body Gaussian form 



Consider the two-body Gaussian potential function in cylindrical coordinates 



The critical first step in the separation method for harmonic-oscillator matrix 
elements is to write the potential itself in a form where the dependence on the 
coordinates ri and r2 has been explicitly separated. We will therefore write 
this two-body function as a sum of one-body functions in the two coordinates. 
Note that the resulting sum will contain and infinite number of terms, while 
the matrix elements of the potential will be limited to a finite sum, thanks to 
properties of the harmonic-oscillator functions. 



B.l Cartesian component 



The radial and Cartesian components of the potential can be expanded inde- 
pendently. We begin with the Cartesian term and postulate 



V {zi,Z2) = e 



choosing for the expansion the functions 



(B.l) 



We will now show that 



(B.2) 



where the coefficients Kz and are given by Eqs. (IB.5I 1 and ( 1B.6I1 . respec- 
tively. 

The exponential function in z'^ in front of the harmonic-oscillator function on 
the left-hand side has been added for computational convenience, as we shall 
see. Then, by orthogonality of the harmonic-oscillator functions, we have 
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OO 



X 



from which we obtain an explicit expression for the weight function {z\ \ hz), 



fn. {zu hz) = / rfz2e"^'/(''^)$„. (Z2; V (Zi, z^) 

J —CO 

Completing the square, we write 



(B.3) 



_f| _ {Zi - Z2f 
hi /X^ 

where we have defined 



hi 



and the integral becomes 



fn. izi;hz)=Afn,exp 



Zl 



GzJ \fi 



X J dz2ex.p — ^ 



fx fx J 



Z2 



(B.4) 



Making the substitutions x = G^J'^Z2lfi, y = G ^^I'^Zxj fx, a = G ^^/"^ fi/hz, the 
remaining integral can be evaluated using Eq. 7.374(8), p. 837 in [8], 



(zi; 6,) = fxG-^I'Mny' (l - exp [- (G. - 1) y= 

ay \ 



xH„ 



After some straightforward algebra and re-grouping of terms, this can be writ- 
ten as 
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or, identifying the term in the square brackets with a harmonic-oscillator func- 
tion with parameter 0^%^ (note the extra factor G^^ needed to get the proper 
normalization constant [py^h-^ ), 



where 



TT/i 



1/2 



G 



(B.5) 
(B.6) 



B.2 Radial component 



For the radial component of the Gaussian potential, we write 



V{pi,p2) = e 



-(pi-p2)2//.2 



OQ CO 



= ^ J2 fnr,A{pUVi;b±)^nr,A{p2,^2;bA 
nr=0 A=— oo 



where we have chosen 




We will then show that 



(B.7) 



(B.8) 



where the coefficients and \2nr+\h\ ^ire given by Eqs. (IB.llll and (IB.12p . 
respectively. 



By orthogonality of the harmonic-oscillator function we then have 



27r '^2 

P2dp2 I ^ dip2e~^e-^P'-P'^"/^"^^^^,Aip2,ip2;b, 
Jo 
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This integral can be evaluated in a straightforward way by transforming to a 
Cartesian coordinate system, and using Eq. (lA.lSp . 



2 I 2 
^2+^2 



^2 I, ,2 



/OO POO -2 

dx2 / dy2e ''^ 

2nr+|A| 

rfy2e"^e-(^-J^^)'/^'<l>„^ (2/2; &±) 



2nr+|A| 



The integrals in the square brackets are precisely those appearing in Eq. (1B.3I1 . 
and they are given by Eq. (1B.2I ) 



2nr+|A| 



nr,A 



nt,=0 



2nr+|A| — riy 



(B.9) 



where 



1 + ^ 



1/2 



and Eq. ()B.9p can be further reduced to 

(Pi, ^1; &±) = i^±A2„.+|A|e-^?/('^-'''-) 

2nr+|A| 

E C:;^^$2n.+|Ahn, (xi; 6^) (z/i; Gf &J 

Jly=0 

Finally, using Eq. (lA.lSP again to return to polar coordinates, we get 



(B.IO) 

(B.ll) 
(B.12) 



/n.,A (Pi, y^i; &±) = i^±A2n.+|A|e'^'/('^^''^)$„,,A (pi, y^i; G 



1/2, 
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C Product of harmonic-oscillator functions 



In this section, we will express the product of two harmonic-oscillator functions 
in terms of a sum of single oscillator functions. These results will be partic- 
ularly useful in evaluating integrals where the integrand includes products of 
harmonic-oscillator functions. 



C.l Product of Cartesian harmonic-oscillator functions 



In this section, we derive the form 



^u, {x- h) {x; h) = , Tl^k2^k {x- h) 

fc=|fci-fc2|,2 



(C.l) 



for the Cartesian harmonic-oscillator functions of Eq. (JT]), with the coefficients 
T^^k^ given by Eq. (ED. 



Using the generating function in Eq. ( lA.ll ). we write for any arbitrary variables 
ti and ^2, 



g-t2+2tia;/fe-a;2/(2fe2)g-t|+2t2a;/fe-a;2/(2fe2) 
CO 2^1/2 



A:i=0 V^l- 



tT^k, {x- h) 



X 



00 2'=2/2 

k2=0 



(C.2) 



With the intent of manipulating the left-hand side of this equation into a form 
similar to the left-hand side of Eq. (lA.lll . we write 



^g-(il+t2)^+2(ti+t2)x/f)-xV(2fe2)^2tit2-a;2/(262) 

Using Eq. ( lA.ll) . this becomes 
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LHS = e2*i*-^V(2.^) ^ E ^ (^1 + $ . (x; h) 



k=0 ^ 
oo 2*^/2 



k=0 



xE 

g=0 



/ 



fe=0 

k f h \ °° OP 

X E E ^tnr^-^ 



5=0 yq J p=o P- 

We can also group the terms in the right-hand side of Eq. (|C.2l) . 



(C.3) 



oo oo 2('^i+'^2)/2 

RHS = bV^Y.Il rrwi ^1^2 ^k. i^: b) $ (x; b) 



(C.4) 



Now we equate powers of ti and ^2 between Eqs. (IC.3I) and ( lC.4p . We find that 
we must make the identifications 



q + p = ki 
k + p- q = k2 

which lead to 



p={ki + k2-k) /2 
q = {ki-k2 + k) /2 



so that Eq. (1C.3P can be written 



ki+k2 



2k/2 



$fc (x; b) 



00 00 

xE E 

fci=0A;2=0 



fe = |fcl-fc2|,2 

k \ 2iki+k2-k)/2 
ki-k2+k 



^ fci+fc2-A: ^ 1 



±kx.k2 



(C.5) 



Note that the limits and step size for the summation over k are dictated by 
the need to keep the arguments of the factorials non-negative. In particular, 
the "2" appearing in the lower limit of the sum over k indicates that the index 
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should be incremented by steps of 2. Direct comparison of Eqs. (1C.4I1 and 
(1C.5I1 now yields 



fci+fc2 



2{'=i+'s2)/2y^ 



|fci-fc2|,2 2 )■ [ 2 )■ [ 2 ) 



-^k (x; b) 



2{ki+k2)/2 

b^^===-<l>k, {x; b) <^k2 {x; b) 



which leads to 



g-x2/(262) ki+k2 

^kAx;b)^k,{x;b)= Tl^k2^k{x-b) 

Ib^ fc=|fci-fc2|,2 



where 



rpk 

^ I ki-k2+k 



1 ^ k2~ki+k ^ \ ^ ki+k2-k ^\ 



(C.6) 



C.2 Product of radial harmonic-oscillator functions 



Here, we obtain the relation 



g-pV(262) '^1,2 

^n^M (P) b) ^n2,k2 {P, b) = ^— J2 Tntkv!^lk2^nM+k2 (P, V\ b) 

V™ n=0 



_(C.7) 

between the harmonic-oscillator functions in polar coordinates defined in Eq. 
(IH) . The expansion coefficients T^f^^n^ fc2 defined by Eq. (1C.9I) . 



Starting from the generating function in Eq. (IA.3I) . and for arbitrary vectors 
ti and t2 



^_p^+2p-u/b-py(2b2)^-ti+2p-t2/b-py{2b2) 

oo oo 



2 51 A/;,,|fci|Xni,fci (ti) (p. b) 

ki=—oo ni=0 

/— oo oo 
,\k2\Xn2,k2 (h) $n2,fc2 (p, b) 



^2=— oo n2=0 

The left-hand side can be written 



(C.8) 
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LHS = e-(*'i+**2)^+2p-(*'i+**2)/fc-pV(262)g2fi-r2-p2/(2f>2) 

Using Eq. HA.Sh again to expand the first exponential, we get 



' ^ fc=-oon=0 

and using Eq. ( IA.37P to absorb the remaining exponential, 

/— oo oo 

k=-oon=0 
oo oo oo oo oo oo 

X E E E E E E 

mi=0 m2=0 pi=0 fci=— oo P2=0 A;2=— oo 



^Xpi+mi,fci (^l) Xp2+ni2,fc2 ( 



Comparing with the right-hand side of Eq. (IC.SP for arbitrary vectors ti and 

— * 

we make the identifications 

Pi + mi = Hi 
P2 + m2 = n2 

and write the left-hand side as 



/— oo oo 

^ fc=-oon=0 
oo oo oo oo oo oo 



X 5^ ^ni-mi,ki,mi;n2-m2,k2,m2 

mi=0 m2=0 711=0 fci=~oo 712=0 fc2=— oo 

xXniM {ti, b) Xn2,k2 {h; b) 
Comparing again with the right-hand side of Eq. flC.81) . we readily deduce 



TT 

2-^ni,|fci|A/;2,|fc2|*ni,fci (P, b) $„,,fc, (p, ^; h) 



/— oo oo 



fc=— OO n=0 

oo oo 



/ y / J n\—m\,k\,m\;n2—m,2,k2,m2 
mi=0 1712=0 

The sum over k disappears because of the Kronecker-delta function inside the 
D coefficient in Eq. (IA.38I) restricting the value of k to ki + k2, and the sum 
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over n is cut off at n = ni,2, because of tlie other Kronecker-delta function in 
Eq. (IA.38I1 restricting its value. Therefore, 
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( 2 J\fn^\ki+k2\ 
^ n=0-^ni,|fci|-^"2,|fc2| 



00 00 



/ y ni— r?ii,A;i,r?ii;n2— rri2ifc2i™2 
mi=0 m2=0 



fcl+fc2 



which we write as 



(P, &) *n2,fe2 (P, = T=T— Tn!',C^lM^^M+k2 (P, 



The coefficients T^^^" ^^^^12 fc2 ^'^^ obtained from Eq. (IA.38I1 , being careful to make 
the substitutions rii ^ rii — mi and n2 — > ^2 — m2 (and therefore, according 
to Eq. (lA.29p . ni 2 — > ni^2 — rrii — m2 as well). Then, 



rpn,ki+k2 
n\,k\;n2,k2 



\n-i+n2-n 




^ n,\n2\{n+\h + h\)\ J^.^W 
I 



ni,2 - nil -7X12 
n 



X 



(ni,2 + \ki + k2 


— mi — 


m2)! 


{ni + 


- mi)! (^2 + /i;2 


- ^2)! 



or, in more compact notation. 



pn,kx+k2 
- n\,k\;n2,k2 



-1) 



ni+n2—n 



n\ (ni + |A;i|)!(n2 + |A:2|)! 
\\ ni!n2! (n+ |/ci + A;2|)! 



ni n2 



r f^n,ki+k2 
^ ^ ""<»ii,2-mi-m2'--'ni,fci,mi;ra2,A:2,ni2 
mi=0 r?i2=0 



(C.9) 



with 
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ni,ki,mi;n2,k2,m2 V / 




(^1,2 + 


\ki + /C2 


— mi — 


1712)1 


{ni + \ki\ 


- mi)! (n2 + \k2 


- 1712)1 



Note again that the Kronecker-delta function 5„<„^ 2-mi-m2 ensures that we 
always have n < ni^2, which we used to limit the sum over n in Eq. (|C.7I) . 



D Formalism for large oscillator shell number 

In this section, we derive the result in [9], 

^6-V2 r(e-ni)r(e-n2)r(e-n) 



(^1 l/nl n2) 



V27r5/2 z^^Jn\nl\n2\ 

X2F1 (-721,-^2; + n + (D.l) 



with ^ given by Eq. (10.60 and 2; by Eq. (ID.lOp . for the numerically accurate 
calculation of the matrix element (ni ^2) in Eq. (fT2l) when large oscillator- 
shell numbers are involved. Note that our result differs slightly from [9] in that 
a "J)-i/2" factors appears in Eq. ( ID. II) instead of "6^/^" (see discussion at the end 
of this section). The formula in Eq. ( ID. II ) is preferred to the one in Eq. ( flQl ) 
for large oscillator-shell numbers, because the latter requires the evaluation 
of a sum of products of large (T) and small (/) coefficients, which can be 
numerically unstable. We also obtain the corresponding matrix elements in 
Eq. jmi) 



^ n, = \n~J^ -n~P\ ,2 



where the coefficients F",^) are defined by Eq. ( ID. 121) . 

n}, ,n\ 

Starting from the definition, 



we use the generating function, Eq. (lA.ll) . to integrate the product of three 



harmonic-oscillator functions with the Gaussian factor. This produces 
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^-tf-tl-t^ / ^2{ti+t2)z/b+2tz/B-iyz^ 



— oo 



oo oo oo 



^ E E E Cn,,n,,ntrtrt'' (^1 \fn\ U^) (D.3) 
ni=0 712=0 n=0 



where 



3 = 0^'% 



1 1 



-'r),i,n2,n 



The integral in the left-hand side of Eq. (ID.Sp is easily evaluated by completing 
the square in the exponential, giving 



LHS=^^e~'"-'^'' +^ (D.4) 
where 

tl + t2 , t 

After some simplification, Eq. (ID. 41) takes the form 
LHS=^exp{[a (ti + ^2) - t]^ C + 2^2} 

with 



G + 1 

which we expand as a series 
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LHS 



\« oo 



1=0 p=0 -f^ 



oo CO 2p 2p— q r)i 



^ ^ ^ ^ 



2* 2p \ 2p- q 



2p-q /•pj.s+i.2p-q-s+i,q 



comparing with the right-hand side of Eq. (ID.3I1 , we make the identifications 



s + i = ni^s = ni — i 
o , • ^ ni+n2 + q . 

q = n 

Note that this implies rii + n2 + n must be even, and the summation over i 
terminates after a finite number of terms, although we will let it run up to oo 
for notational convenience, letting the factorial terms implicitly truncate the 
sum. Then we have 



Cm, 112, n {ni \fn\ n2) 



— |'_q;)"i+"2 ^(ni+n2+n)/2 



ni + n2 + n — 2i 
n 



oo 

i=0 
7f 



^ rii + n2 — 
^ rii-i y 



i+n2+» 



V 

oo 



f Biz 



[rii + n2 + n — 2i)\ 



fo n\ (ni - i)\ (712 - iV- ( "^+2'^" - ^)-^- V"^^' 
Next, we simplify the ratio of factorials 



(D.5) 



(2p)! _ (ni + n2 + n - 2i)\ 

p\ ^ ni+?i2+ra ^^1 

using the doubling formula for the Gamma function (Eq. 8.335(1) in f^), 



pi VTT V 2 



For convenience, we define 
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^^ „. ^^^^ 

which is a half-integer since we have already noted that ni + n2 + n is even. 
Then p = ^ — i — 1/2^ and 



In order to simplify this further, we derive the following useful identity 

r(i-e + = (^-Or(^-0 

= (^-0 (^-e-i)---(i-Or(i-0 

= (-ir(e-i)---(e-(^-i))(e-Or(i-o 

Similarly, we can write 

r(O = (^-i)---(^-(.-i))(e-0r(e-0 

Therefore, 



and, equivalently. 



r(e-^) = (-ir^^^|^|^ (D.9) 

Thus, Eq. ( ID .Ti l becomes 



(2p)! 2^^ , r(or(i-o 



' (1-0. 

where we have used the Pochhammer symbol 



r (x + n) 

= = X (x + 1) ■ ■ ■ (x + n - 1) 

Returning to Eq. (ID. 51) . we replace the [rii — i)\ and {n2 — i)\ terms with 
Pochhammer symbols as well using Eq. ( ID. 91 ) with ,^ ^ ni + 1 to write 
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{rii — i)\ = T {rii — i + 1) 

r(ni + i)r(-ni) 



i-iy 
{-ly 



T{-ni+i) 
rill 



and similarly for (^2 — Then, Eq. (ID.SP yields 



which we express as a hypergeometric function, as defined in [8] Eq. 9.100 (see 
also section 9.14(2) in [8] for the notation in terms of a generalized hypergeo- 
metric function), 



2g-VVr (0 i-ar^-^ ^(ni+n.+.)/2 

|/n| ?^2) = ^ — ^ ■ — 2F1 (-ni, -ria; 1 - 4; ^) 



where 



^ 1 + ^ (D-10) 



2a\ 262 
Simplifying further, we find 

{n, \U\ n,) = (_l)(-+--")/2 / ^^^] 

X2Fi{-ni,-n2]l-i]z) (D.ll) 

Comparing with Eq. (3) in [9], we note that the hypergeometric function is 
evaluated at 1 — 2; rather than z in that paper. In order to make a direct 
comparison with [9], we use Eq. 9.131(2) in [8], 



.Fi (-ni, -n2] ~ni - n2 + z) 



2F1 {-rii, -n2; I - C.; z) 

r(i-or(i-e + ni + n2) 
r(i-e + ni)r(i-e + n2) 

r(-ni)r(-n2) 

X 2F1 (1 - ^ + m, 1 - C + n2; 2 - ^ + m + n2; 1 - 2:) 



l-5+ni+n2 r(l-Or(-ni- 712+^-1) 
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The second term vanishes because of the Gamma functions with negative- 
integer (or zero) arguments in the denominator. We can simplify the third 
argument of the hypergeometric function in the first term to 



-rii — n2 + n + I 



+ n + 1 



Thus, 



2F1 (-ni, -ria; I - ^; 



r (1 - e + Til) r (1 - ^ + 722) 
X 2F1 {-rii, -n2; -^ + n + l;l- z) 



Next, we use Eq. ()D.8[1 to re- write the Gamma functions in the denominator, 

2F, (-n„ -.2, i-^;z) = (-1) r(i-or(Or(0 

X 2F1 (-ni, -^2; + n + 1;1 - z) 
Substituting this expression into Eq. f ID.llll gives 

X 2-Fi (-^1, -n2; -^ + n + l;l - z) 
Finally, we use Eq. 8.334(3) in {8| to write 



{n, |/„| n2) = (-!)(-— )/^ / r(e-n0r(e-n2)r(^-n) 

^J2hy/T{ ^/n\ni\n2\T^ [-ly 
X 2-Fi (-JT-i, -712; + n + 1; 1 - 2;) 

^ r (^ - m) r (e - ^2) r (e - n) 

X 2-Fi (-721, -n2; n 1; 1 - 2;) 

This result is nearly identical to Eq. (3) in [9], after properly adjusting for 
the choice of variable names, the only minor difference being the oscillator 
parameter which appears as in the present work, and 6^^^ in |9]. However, 
dimensional analysis favors the form, as the matrix element {rii \ fn\n2) 
must carry dimensions of length to the 1/2 power, according to its definition in 
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Eq. (fT2ll . In closing, we use Eq. (ID.ip to write the expression for the two-body 
matrix element (corresponding to Eq. (fTOl) in the large oscillator-shell limit), 



1/(2) _ ^ \ ^ rpriz rnriz 



where 



(i) (fc) 



X 



Fi (-n«, -ni'^); + + 1; 1 - ^) (D.12) 



E Angular integral 

We wish to evaluate the radial part of the matrix-element integral 

, ^ roo /'27r roo r2n 

= Pidpi dipi p2dp2 dip2 

^*n('\A« (^1' '^n'^'-AO) ^^2, ¥^2; &±) 

Xe-(^^i-^^^)VM^$^,,)^^^^^ (Pl, <^i; &±) '^'„(0_^(o (P2, v^2; &±) 

numerically, where the harmonic-oscillator functions are defined in Eq. dH). 
By rotational invariance of the Gaussian potential, we have 

_A» - A(^') + A^'^) + A(') = 

The angular integrals over (pi and (p2 are particularly problematic because of 
their oscillatory nature. Therefore, we focus on those integrals and introduce 
the function 

1 /'27r /'27r 
ek{x) = 2/ d^i ^^^^ik{^,~^,)^2xcos(^,-^,) 

(27r) Jo Jo 
so that we may write 
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/ N I OO coo 
(r) 



Pidpi P2dp2e |^«| (pi;6±)$„u)_|^0)| (P2;&±) 

*np\|AW| (/^l! ^J-) ^n«,|A(0| (Ps; 6_AW+A(fc) ^^^^ (^-2) 

We simplify Eq. (lE.ip using the generating function for the Bessel function, 
given in Eq. flA.611 . with z = —2ix and (p = (pi — {p2, 

oo 
n=— oo 

from which the integral in Eq. (IE. 10 yields 

e,.(a;)=^l'lj|fe|(-2zx) 

From the series expansion of the modified Bessel function of the first kind, Eq. 
8.445 in [8], we get 

e, (x) = (-i)i'i i\k\{-2x) 

oo 



We find that the series in Eq. (|E.3l) is extremely well converged if we include 
terms up to m such that 



m\ (m + \k\)\ 



< e 



where e = IQ-'^'^o-N^u^d/s ^qj, calculation in up to Nq oscillator shells and 
-^quad quadrature points. The remaining integrals over pi and p2 in Eq. f lE.2l ) 
were evaluated by Gauss-Laguerre quadrature. 
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